﻿///////////////////////////////////////////////////////////////////////////////
//
//  This file is part of MathLib.NET.
//
//  This library is free software; you can redistribute it and/or
//  modify it under the terms of the GNU Lesser General Public
//  License as published by the Free Software Foundation; either
//  version 2.1 of the License, or (at your option) any later version.
//  
//  This library is distributed in the hope that it will be useful,
//  but WITHOUT ANY WARRANTY; without even the implied warranty of
//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
//  Lesser General Public License for more details.
//  
//  You should have received a copy of the GNU Lesser General Public
//  License along with this library;  If not, see 
//  <http://www.gnu.org/licenses/>.
//
///////////////////////////////////////////////////////////////////////////////

using System;
using ILNumerics;
using ILNumerics.BuiltInFunctions;

// Complex functions
// =================
//
// abs - Absolute value and complex magnitude 
// angle - Phase angle 
// complex - Construct complex data from real and imaginary components 
// conj - Complex conjugate 
// cplxpair - Sort complex numbers into complex conjugate pairs 
// i - Imaginary unit 
// imag - Imaginary part of complex number 
// isreal - Check if input is real array 
// j - Imaginary unit 
// real - Real part of complex number 
// sign - Signum function 
// unwrap - Correct phase angles to produce smoother phase plots

namespace MathLib
{
    /// <summary>
    /// Complex functions not provided by ILNumerics.
    /// </summary>
    public partial class MLMath
    {
        //
        // --- Absolute value and complex magnitude ---
        //

        /// <summary>
        /// Absolute value and complex magnitude.
        /// </summary>
        /// <param name="arg">Argument to return absolute value of</param>
        /// <returns>Absolute value of argument</returns>
        public static ILArray<double> abs(ILArray<double> arg)
        {
            return ILMath.abs(arg);
        }

        /// <summary>
        /// Absolute value and complex magnitude.
        /// </summary>
        /// <param name="arg">Argument to return absolute value of</param>
        /// <returns>Absolute value of argument</returns>
        public static ILArray<float> abs(ILArray<float> arg)
        {
            return ILMath.abs(arg);
        }

        /// <summary>
        /// Absolute value and complex magnitude.
        /// </summary>
        /// <param name="arg">Argument to return absolute value of</param>
        /// <returns>Absolute value of argument</returns>
        public static ILArray<int> abs(ILArray<int> arg)
        {
            return ILMath.abs(arg);
        }

        /// <summary>
        /// Absolute value and complex magnitude.
        /// </summary>
        /// <param name="arg">Argument to return absolute value of</param>
        /// <returns>Absolute value of argument</returns>
        public static ILArray<double> abs(ILArray<complex> arg)
        {
            return ILMath.abs(arg);
        }

        /// <summary>
        /// Absolute value and complex magnitude.
        /// </summary>
        /// <param name="arg">Argument to return absolute value of</param>
        /// <returns>Absolute value of argument</returns>
        public static ILArray<float> abs(ILArray<fcomplex> arg)
        {
            return ILMath.abs(arg);
        }

        /// <summary>
        /// Absolute value and complex magnitude.
        /// </summary>
        /// <param name="arg">Argument to return absolute value of</param>
        /// <returns>Absolute value of argument</returns>
        public static ILArray<byte> abs(ILArray<byte> arg)
        {
            return ILMath.abs(arg);
        }

        /// <summary>
        /// Absolute value and complex magnitude.
        /// </summary>
        /// <param name="arg">Argument to return absolute value of</param>
        /// <returns>Absolute value of argument</returns>
        public static ILArray<char> abs(ILArray<char> arg)
        {
            return ILMath.abs(arg);
        }

        /// <summary>
        /// Absolute value and complex magnitude.
        /// </summary>
        /// <param name="arg">Argument to return absolute value of</param>
        /// <returns>Absolute value of argument</returns>
        public static ILArray<long> abs(ILArray<long> arg)
        {
            return ILMath.abs(arg);
        }

        /// <summary>
        /// Absolute value and complex magnitude.
        /// </summary>
        /// <param name="arg">Argument to return absolute value of</param>
        /// <returns>Absolute value of argument</returns>
        public static ILArray<short> abs(ILArray<short> arg)
        {
            return ILMath.abs(arg);
        }


        //
        // --- Complex conjugate ---
        //

        /// <summary>
        /// Complex conjugate.
        /// </summary>
        /// <param name="arg">Argument to return conjugate of</param>
        /// <returns>Conjugate of argument</returns>
        public static ILArray<complex> conj(ILArray<double> arg)
        {
            return ILMath.conj(ILMath.tocomplex(arg));
        }

        /// <summary>
        /// Complex conjugate.
        /// </summary>
        /// <param name="arg">Argument to return conjugate of</param>
        /// <returns>Conjugate of argument</returns>
        public static ILArray<fcomplex> conj(ILArray<float> arg)
        {
            return ILMath.conj(ILMath.tofcomplex(arg));
        }

        /// <summary>
        /// Complex conjugate.
        /// </summary>
        /// <param name="arg">Argument to return conjugate of</param>
        /// <returns>Conjugate of argument</returns>
        public static ILArray<complex> conj(ILArray<int> arg)
        {
            return ILMath.conj(ILMath.tocomplex(arg));
        }

        /// <summary>
        /// Complex conjugate.
        /// </summary>
        /// <param name="arg">Argument to return conjugate of</param>
        /// <returns>Conjugate of argument</returns>
        public static ILArray<complex> conj(ILArray<complex> arg)
        {
            return ILMath.conj(arg);
        }

        /// <summary>
        /// Complex conjugate.
        /// </summary>
        /// <param name="arg">Argument to return conjugate of</param>
        /// <returns>Conjugate of argument</returns>
        public static ILArray<fcomplex> conj(ILArray<fcomplex> arg)
        {
            return ILMath.conj(arg);
        }

        /// <summary>
        /// Complex conjugate.
        /// </summary>
        /// <param name="arg">Argument to return conjugate of</param>
        /// <returns>Conjugate of argument</returns>
        public static ILArray<complex> conj(ILArray<byte> arg)
        {
            return ILMath.conj(ILMath.tocomplex(arg));
        }

        /// <summary>
        /// Complex conjugate.
        /// </summary>
        /// <param name="arg">Argument to return conjugate of</param>
        /// <returns>Conjugate of argument</returns>
        public static ILArray<complex> conj(ILArray<char> arg)
        {
            return ILMath.conj(ILMath.tocomplex(arg));
        }

        /// <summary>
        /// Complex conjugate.
        /// </summary>
        /// <param name="arg">Argument to return conjugate of</param>
        /// <returns>Conjugate of argument</returns>
        public static ILArray<complex> conj(ILArray<long> arg)
        {
            return ILMath.conj(ILMath.tocomplex(arg));
        }

        /// <summary>
        /// Complex conjugate.
        /// </summary>
        /// <param name="arg">Argument to return conjugate of</param>
        /// <returns>Conjugate of argument</returns>
        public static ILArray<complex> conj(ILArray<short> arg)
        {
            return ILMath.conj(ILMath.tocomplex(arg));
        }


        //
        // --- Angle in complex plane ---
        //

        /// <summary>
        /// Angle in complex plane.
        /// </summary>
        /// <param name="arg">Argument to return absolute value of</param>
        /// <returns>Absolute value of argument</returns>
        public static ILArray<double> angle(ILArray<double> arg)
        {
            return ILMath.zeros(arg.Dimensions);
        }

        /// <summary>
        /// Absolute value and complex magnitude.
        /// </summary>
        /// <param name="arg">Argument to return absolute value of</param>
        /// <returns>Absolute value of argument</returns>
        public static ILArray<float> angle(ILArray<float> arg)
        {
            return ILMath.tosingle(ILMath.zeros(arg.Dimensions));
        }

        /// <summary>
        /// Absolute value and complex magnitude.
        /// </summary>
        /// <param name="arg">Argument to return absolute value of</param>
        /// <returns>Absolute value of argument</returns>
        public static ILArray<double> angle(ILArray<int> arg)
        {
            return ILMath.zeros(arg.Dimensions);
        }

        /// <summary>
        /// Absolute value and complex magnitude.
        /// </summary>
        /// <param name="arg">Argument to return absolute value of</param>
        /// <returns>Absolute value of argument</returns>
        public static ILArray<double> angle(ILArray<complex> arg)
        {
            return ILMath.atan2(ILMath.real(arg), ILMath.imag(arg));
        }

        /// <summary>
        /// Absolute value and complex magnitude.
        /// </summary>
        /// <param name="arg">Argument to return absolute value of</param>
        /// <returns>Absolute value of argument</returns>
        public static ILArray<float> angle(ILArray<fcomplex> arg)
        {
            return ILMath.tosingle(ILMath.atan2(ILMath.todouble(ILMath.real(arg)), ILMath.todouble(ILMath.imag(arg))));
        }

        /// <summary>
        /// Absolute value and complex magnitude.
        /// </summary>
        /// <param name="arg">Argument to return absolute value of</param>
        /// <returns>Absolute value of argument</returns>
        public static ILArray<double> angle(ILArray<byte> arg)
        {
            return ILMath.zeros(arg.Dimensions);
        }

        /// <summary>
        /// Absolute value and complex magnitude.
        /// </summary>
        /// <param name="arg">Argument to return absolute value of</param>
        /// <returns>Absolute value of argument</returns>
        public static ILArray<double> angle(ILArray<char> arg)
        {
            return ILMath.zeros(arg.Dimensions);
        }

        /// <summary>
        /// Absolute value and complex magnitude.
        /// </summary>
        /// <param name="arg">Argument to return absolute value of</param>
        /// <returns>Absolute value of argument</returns>
        public static ILArray<double> angle(ILArray<long> arg)
        {
            return ILMath.zeros(arg.Dimensions);
        }

        /// <summary>
        /// Absolute value and complex magnitude.
        /// </summary>
        /// <param name="arg">Argument to return absolute value of</param>
        /// <returns>Absolute value of argument</returns>
        public static ILArray<double> angle(ILArray<short> arg)
        {
            return ILMath.zeros(arg.Dimensions);
        }


        //
        // --- Determine if real (not complex) ---
        //

        /// <summary>
        /// Determines if the supplied argument is a real entity.
        /// </summary>
        /// <param name="obj">Object to test</param>
        /// <returns>true if real, otherwise false</returns>
        public static bool isreal(object obj)
        {
            if (obj is int ||
                obj is uint ||
                obj is char ||
                obj is byte ||
                obj is float ||
                obj is double ||
                obj is ILArray<int> ||
                obj is ILArray<uint> ||
                obj is ILArray<char> ||
                obj is ILArray<byte> ||
                obj is ILArray<float> ||
                obj is ILArray<double>)
                return true;

            return false;
        }


        //
        // --- Correct phase angles ---
        //

        /// <summary>
        /// Correct phase angles to produce smoother phase plots.
        /// </summary>
        /// <param name="angle">Matrix containing phase angles</param>
        /// <param name="tol">Maximum jump tolerance</param>
        /// <param name="dim">Dimension to operate along</param>
        /// <returns>Array containing corrected phase angles</returns>
        private static ILArray<T> unwrap<T, C>(ILArray<T> angle, object tol = null, int dim = 1)
            where T : new()
            where C : INumber<T>, new()
        {
            ILArray<T> unwrapArray = new ILArray<T>(angle.Dimensions.ToIntArray());

            // Create arithmetic calculator
            C calc = new C();

            // Determine tolerance to use
            T tolerance;
            if ((tol is uint || tol is int || tol is float || tol is double) && (double)tol > Math.PI)
                tolerance = calc.FromDouble((double)tol);
            else
                tolerance = calc.FromDouble(Math.PI);

            T twoPi = calc.FromDouble(2.0 * Math.PI);

            // Calculate strides and counts
            int numTransforms, majorTransformStride, numPoints, pointStride;
            ArrayHelper.CalcArrayAccessParams<T>(angle, ref dim, out numTransforms, out majorTransformStride, out numPoints, out pointStride);

            for (int maj = 0; maj < numTransforms; maj++)
            {
                for (int min = 0; min < pointStride; min++)
                {
                    int startIdx = maj * majorTransformStride + min;

                    T lastAngle = unwrapArray.InternalArray4Experts[startIdx];
                    T offset = calc.FromInt(0);
                    for (int i = 0; i < numPoints; i++)
                    {
                        T a = angle.InternalArray4Experts[startIdx + i * pointStride];
                        if (calc.Le(calc.Abs(calc.Sub(a, lastAngle)), tolerance))
                        {
                            // Within tolerance
                            unwrapArray.InternalArray4Experts[startIdx + i * pointStride] = calc.Add(a, offset);
                        }
                        else if (calc.Gt(a, lastAngle))
                        {
                            // Outside tolerance and greater than last
                            while (calc.Gt(calc.Add(a, calc.Sub(offset, lastAngle)), tolerance)) offset = calc.Sub(offset, twoPi);
                            unwrapArray.InternalArray4Experts[startIdx + i * pointStride] = calc.Add(a, offset);
                        }
                        else
                        {
                            // Outside tolerance and less than last
                            while (calc.Lt(calc.Add(a, calc.Sub(offset, lastAngle)), tolerance)) offset = calc.Add(offset, twoPi);
                            unwrapArray.InternalArray4Experts[startIdx + i * pointStride] = calc.Add(a, offset);
                        }
                        lastAngle = a;
                    }
                }
            }

            return unwrapArray;
        }

        /// <summary>
        /// Correct phase angles to produce smoother phase plots.
        /// </summary>
        /// <param name="angle">Matrix containing phase angles</param>
        /// <param name="tol">Maximum jump tolerance</param>
        /// <param name="dim">Dimension to operate along</param>
        /// <returns>Array containing corrected phase angles</returns>
        public static ILArray<double> unwrap(ILArray<double> angle, object tol = null, int dim = 1)
        {
            return unwrap<double, MLDouble>(angle, tol, dim);
        }

        /// <summary>
        /// Correct phase angles to produce smoother phase plots.
        /// </summary>
        /// <param name="angle">Matrix containing phase angles</param>
        /// <param name="tol">Maximum jump tolerance</param>
        /// <param name="dim">Dimension to operate along</param>
        /// <returns>Array containing corrected phase angles</returns>
        public static ILArray<complex> unwrap(ILArray<complex> angle, object tol = null, int dim = 1)
        {
            return unwrap<complex, MLComplex>(angle, tol, dim);
        }

        /// <summary>
        /// Correct phase angles to produce smoother phase plots.
        /// </summary>
        /// <param name="angle">Matrix containing phase angles</param>
        /// <param name="tol">Maximum jump tolerance</param>
        /// <param name="dim">Dimension to operate along</param>
        /// <returns>Array containing corrected phase angles</returns>
        public static ILArray<float> unwrap(ILArray<float> angle, object tol = null, int dim = 1)
        {
            return unwrap<float, MLFloat>(angle, tol, dim);
        }

        /// <summary>
        /// Correct phase angles to produce smoother phase plots.
        /// </summary>
        /// <param name="angle">Matrix containing phase angles</param>
        /// <param name="tol">Maximum jump tolerance</param>
        /// <param name="dim">Dimension to operate along</param>
        /// <returns>Array containing corrected phase angles</returns>
        public static ILArray<fcomplex> unwrap(ILArray<fcomplex> angle, object tol = null, int dim = 1)
        {
            return unwrap<fcomplex, MLFComplex>(angle, tol, dim);
        }
    }
}
